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Abstract 



Questions surrounding the spatial disposition of particles in various condensed-matter 
systems continue to pose many theoretical challenges. This paper explores the geometric 
availability of amorphous many-particle configurations that conform to a given pair correlation 
function g{r). Such a study is required to observe the basic constraints of non-negativity for 
g{r) as well as for its structure factor S{k). The hard sphere case receives special attention, 
to help identify what qualitative features play significant roles in determining upper limits 
to maximum amorphous packing densities. For that purpose, a five-parameter test family 
of g^s has been considered, which incorporates the known features of core exclusion, contact 
pairs, and damped oscillatory short-range order beyond contact. Numerical optimization over 
this five-parameter set produces a maximum-packing value for the fraction of covered volume, 
and about 5.8 for the mean contact number, both of which are within the range of previous 
experimental and simulational packing results. However, the corresponding maximum-density 
g{r) and S{k) display some unexpected characteristics. These include absence of any pairs at 
about 1.4 times the sphere collision diameter, and a surprisingly large magnitude for S{k = 0), 
the measure of macroscopic-distance-scale density variations. On the basis of these results, 
we conclude that restoration of more subtle features to the test-function family of g^s (i.e. 
a split second peak, and a jump discontinuity at twice the collision diameter) will remove 
these unusual characteristics, while presumably increasing the maximum density slightly. A 
byproduct of our investigation is a lower bound on the maximum density for random sphere 
packings in d dimensions, which is sharper than a well-known lower bound for regular lattice 
packings for d > 3. 



1 Introduction 



Over a broad range of length scales, many-particle systems exhibit a rich variety of structures 
with varying degrees of long-range order, spanning from crystals, quasicrystals, and polycrys- 
tals to amorphous solids and liquids. Consequently, it is natural to focus attention on the 
statistical mechanics of the arrangement of the particles. In the case of a macroscopic system 
containing a large number of particles, a full configurational description of that system 
usually is neither feasible, desirable, nor necessary. For most practical purposes, it suffices to 
determine, or to describe, the distribution functions of low orders n <^ N. Conventionally, this 
information is conveyed in the form of correlation functions. For statistically homogeneous 
systems consisting of identical spherical particles in a volume V, these correlation functions 
are defined so that p^g^'^\vi^ r2, . . . , r„) is proportional to the probability density for simulta- 
neously finding n particles at locations ri, r2, . . . , r„ within the system,0 where p = N/V is the 
number density. With this convention, each g^'^^ approaches unity when all particle positions 
become widely separated within V. 

The present study concerns the special circumstances for which the constituent parti- 
cles are spherically symmetric and identical, and the system is statistically homogeneous and 
isotropic. These conditions can be satisfied if the system contains a single fiuid or amorphous 
solid phase. The correlation function of primary interest is g^'^^r^), depending configura- 
tionally just on scalar pair distance ru, and thus specifying how many pair distances of a 
given length occur statistically within the system. The third-order function g^^\ri2,ri2,r23) 
reveals how these pair separations are linked into triangles. This additional information strictly 
speaking cannot be inferred from the knowledge oi g = g^'^^ alone, although the Kirkwood su- 
perposition approximationlili presumes to fill that knowledge gap. The fourth-order function 
(yf^^'' controls the assembly of triangles into tetrahedra, and is the lowest-order correlation 
function that is sensitive to chirality of the medium. 

On account of their probability interpretation, all of the ^f*^"^ must be non-negative; in 
particular, for all r > we must have 

g{r) > 0. (1) 
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In addition to this fundamental constraint, g{r) is also subject to another basic inequality that 
arises from its connection to density fluctuations. This concerns the behavior of the structure 
factor defined thus:! 

S{k) = 1 + p J exp(— ik ■ r)[g{r) — l]dr 

/■°° r sin fcr . , , ^ , , _ , 

= l + inp -^—[g{r) - l]dr. (2) 

The second line assumes that we are treating three-dimensional systems. The second funda- 
mental constraint is the non- negativity of S{k), i.e., 

S{k) > 0, (3) 

which must be obeyed for all real values of k. It should be noted that (|ip and are not at 
all restricted to states of thermal equilibrium, but are more general. It is currently unknown 
if these necessary conditions (|I]) and @ are also sufficient to guarantee that any function 
satisfying them is actually the pair correlation function for a realizable many-particle system; 
however, no counterexamples are currently known. 

Two recent studied'! have examined the theoretical possibility of controlling pair inter- 
actions in an isothermal many-particle system over a non- vanishing density range, starting at 
p = 0, in such a way that g{r) remains invariant over that range. The cases examined have 
assigned forms to the invariant g that were the zero-density limits appropriate for rigid rods, 
disks, and spheres,! and for the hard core plus square well pair potential.! In each of these 
examples an upper terminal density p* could be identified such that over 

< p < p* (4) 

the g invariance could indeed be maintained. However, crossing p* would cause the S{k) 
inequality (H) to be violated at some k. 

The principal objective of the present project is to apply the (7-invariance technique to the 
still-challenging problem of random packings of spheres. It has now been well established that 
the old concepts of "random loose packing" and "random close packing" are ill-defined.0 In- 
stead, a non-trivial density range exists over which irregular packings of various types (locally, 



collectively, or strictly jammedEJ) can be formed, the preferred densities and local structures 




of which depend on the preparation algorithm. Our objective has been to study the logical 
connections between qualifying rigid-sphere g^s and the maximum corresponding densities, 
emerging as the terminal p*'s. 

The following Section II introduces a parametric family of qualitatively reasonable, but 
functionally elementary, pair correlation functions for the amorphous-state sphere spatial ar- 
rangement problem. This family contains a set of five adjustable parameters, whose values 
must of course be consistent with the two basic inequalities (|I|) and (^). Section III describes 
a numerical search procedure over these parameters, and its results, the goal of which was to 
produce the largest p*. We believe it is significant that even with such a simple parametric 
family of pair correlation functions, p* can come close to that obtained in many experimental 
and simulational preparations of random sphere packings. Section IV offers some interpretive 
remarks stimulated by the numerical results in Section III, and indicates the natural and 
useful directions for future investigation. Finally, in an appendix, we derive a lower bound 
on the maximum density p* for random sphere packings in d dimensions and show that it is 
sharper than a well-known lower bound for regular lattice packings for d > 3. 



Our interest in this paper is to study models in which long-range order is suppressed and 
short-range order is controlled. Information that is available from previous determinations 
of local order in amorphous sphere packings provides useful guidance in choosing a model 
family of functions for the present investigation. In particular, we note that a survey of sev- 



procedures, appear to agree on the presence of some qualitative features. (Figure |l| shows the 
pair correlation function for a dense random packing of spheres as generated by us employ- 
ing the Lubachevsky-Stillinger "compression" protocol.0) Using the sphere-pair distance of 
closest approach as the natural length unit, these g{r) attributes are the following: 



2 Model Family 



eral experimental and computer-simulation 




using distinct packing preparation 
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(i) Obviously, g{r) must vanish for all < r < 1. 

(ii) On account of the jamming, virtually all spheres (a few "rattler" spheres can be present 
as exceptions) are rigidly in contact with neighbors. The number of such contacts must 
average at least 4 to meet the definition of "local" jamming .0 

(iii) For r > 1, g{r) displays finite-amplitude oscillations about unity, that decay to zero 
with increasing r. The length scale of these oscillations is roughly comparable to the 
sphere diameter. 

(iv) A pair of distinctive g{r) peaks appear at distances approximately equal to -\/3 and 
a/4. These are often termed a "split second peak", and appear in modified form for 
amorphous deposits of soft-sphere and attracting particles.0 

(v) As r increases through r = 2, g{r) experiences a discontinuous drop in magnitude. 

It is not the objective of the present work to try to include all of these features slavishly. 
Instead, we have chosen as a first step to represent only attributes (i), (ii), and (iii) by a simple 
parametric function family, and to see how close the largest corresponding density would come 
to the approximate experimental rang e§'0 of "random" packing densities:lli 

1.18 <p< 1.26 (5) 
The equivalent approximate range of covering or packing fractions = 7rp/6 is: 

0.62 < < 0.66 (6) 
Computer-simulation determinations of random packing fractionJ^ lie in the wider range 

0.60 < < 0.68, (7) 

illustrating the fact that the packing densities are protocol-dependent. I' I After the fact, this 
approach should determine how important the remaining attributes (iv) and (v) are for the 
random sphere packing problem. 
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Consequently, we have elected to write g{r) as a linear combination of three portions, 
corresponding respectively to (i), (ii), and (iii) above: 



g{r) = gi{r) + gn{r) + giii{r] 



(8) 



The first involves just the unit step function U: 



gi{r)^Uir-l), 



(9) 



while the second represents the sphere contacts: 

'7 

gu{r) 



S{r-1). 



(10) 



Here, Z is the first of our adjustable parameters, equal to the mean number of contacts 
(coordination number) experienced by each sphere. The third portion contains four more 
adjustable dimensionless parameters {A, B, C, D): 



giii{r) = - exp(-Sr) sin(Cr + D)U{r - 1), 
r 

and is intended to represent approximately the damped oscillation beyond contact. 
The Fourier transforms required to evaluate the structure factor, 

S{k) = 1 + p[Gi{k) + Gn{k) + Gni{k)] 
are expressible entirely in terms of elementary functions: 

Gi{k) — y^exp(— ik • v)[gi{r) — l\dv 

— — f/c cos k — sin k] 
k"^ 

Gn{k) = j exp(-ik • r)[gu{r) - Ijdr 

Z sink 
p k 



(11) 



(12) 



(13) 



(14) 



Gni{k) = j exp(-ik ■ Y)[giii{r) - l]rfr 



27rexp(-S) 
k 



Bcos{k -C-D)-{k-C) siTL{k -C-D) 
B^ + {k-CY 



Bcos{k + C + D) -{k + C) sin(A; + C + D) 

B^ + {k + cy 
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(15) 



The five adjustable parameters Z, A, B, C, and D are subject to some obvious constraints. 
Clearly, we must demand that 

Z >0. (16) 

In addition, exponential increase with r is not permissible, so 

B>0. (17) 

Furthermore, D only needs to span a single period of the trigonometrical factor in which it 
occurs: 

< D < 27r. (18) 

The remaining parameter pair A, C is not entirely free, of course, but must be consistent 
with both inequalities (|l]) and (|^). Our central objective is to search over the five- dimensional 
domain defined by (||), (H), and (p^)-([T8|), for the maximum terminal density p*{Z, A, B, C, D) 
or terminal covering volume fraction 0*(Z, A, 5, C, Z)) that they permit. Under the working 
assumption that this maximum is attained at the boundary of the five-dimensional domain, 
it becomes important to know which constraint or constraints are at issue there, and why. 



3 Numerical Search Procedure and Results 

The problem of finding the maximal packing fraction (f)*{Z, A, B,C, D) can be posed as an 
optimization problem. Specifically, this can be posed as a two-level "min-max" problem: one 
wants to maximize (piZ, A, B, C, D) over the parameters Z, A, B, C, and D with the restrictions 
(p!^)-([T8|) such that the minimum of S{k] Z, A, B, C, D) in the variable k and the minimum of 
g{r; Z, A, B, C, D) in the variable r are both non-negative, i.e., 

max S (19) 

Z,A,B,C,D ^ ' 

such that 

min5(A;;Z,A,5,<^,^) > 0, min c/(r; Z, A, 5, ^) > 0. (20) 

k 

The interval-arithmetic paradig is a global optimization methodology that in principle 
should enable one to obtain exact narrow interval bounds on the maximal packing fraction in 



a computationally efficient manner. We attempted this calculation using the GlobSol Fortran 
90 global optimization libraryS but could not obtain an exact interval solution. This program 
is best suited for finding conventional extrema of simple differentiable functions with simple 
constraints. Our problem is considerably more complex and so we instead used GlobSol to 
find the minimum 

mmS{k;Z,A,B,C,D) 

k 

and then employed a brute-force grid search over Z, A, B, C, and D subject to the aforemen- 
tioned conditions in order to maximize 0. 

The search procedure is implemented for six different cases summarized in Table 1. In the 
first case (Case I), no restrictions on the functions, other than the ones described above, are 
imposed. In the remaining cases, we impose additional restrictions. In particular, it is known 
that dense random packings are typically spatially uniform, implying that S{k = 0) ~ 0. 
Therefore, in some instances, we carry out the search subject to the condition that S{k = 
0) = 0. 

Case I 

Not surprisingly, the least restrictive case yields the largest value of the maximal packing 
fraction: (p* = 0.627. The corresponding values of the parameters are listed in Table 2. Figure 
^ shows the structure factor and pair correlation function. These functions reveal structural 
features that are not characteristic of typical dense random packings obtained either exper- 
imentally or computationally. For example, the structure factor at A; = is unusually high, 
implying significant density fiuctuations in the infinite- wavelength limit. Moreover, at three 
different finite wavelengths (A; ^ 3.1, ~ 6.3, and k ^ 10) the structure factor is essentially 
zero or nearly zero, implying vanishingly small density fiuctuations at these wavelengths. 
Atypically, the pair correlation function g{r) exhibits appreciable oscillations before attaining 
its long-range value at about r = 8. This feature could be the result of polycrystallinity, but 
the fact that g{r) vanishes at r ^ 1.4 (another anomalous characteristic) evidently eliminates 
the possibility that such putative crystallites are face-centered cubic arrangements. Inter- 
estingly, the coordination number Z ^ 5.8 is approximately equal to the value observed in 
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experiments and simulations of typical dense random packings.0a Theoretical arguments have 
been put forth predicting that Z = Q for random packings of identical frictionless spheres in 
three dimensions. 



Case II 



In the second case, we conduct the search under the condition that S{k = 0) = 0. For 
this condition to hold, Z must be given by 



^0-1 



240 



B 



sm{C + D) + (C 



2BC 



cos(C + D) 



■ (21) 



B^ + cy ' ' V B^ + C^ 

This condition also implies that the term in S{k) of order k"^ is zero, and therefore the first 
non-zero term is of order k'^. Both the maximal packing fraction and coordination number 
drop from the first case to the values 0* = 0.46 and Z = 2.3964 (see Tables 1 and 2). Figure 
^ shows the structure factor and pair correlation function. Note the atypical curvature of the 
function g{r) near the contact value. 



Case III 

In the third case, we suppress the damped-oscillating component [giii{r) = Guiik) = 0]. 
Here we find that 0* = 0.41 and Z = 3.1504. 



Case IV 

In the fourth case, we suppress the damped-oscillating component [gin{r) = Gin{k) = 0] 
and we also impose the condition S{k = 0) = 0. This problem can be solved exactly. Here we 
find 




In the Appendix, we obtain the (i-dimensional generalization of this result and show how 
it leads to a lower bound on 0* for the general case in which the third component is not 
suppressed [giii{r) ^ 0, Gni{k) ^ 0]. 



Case V 
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In the fifth case, we suppress the delta-function component {Z = 0). Here we obtain 
(j)* = 0.375 (see Tables 1 and 2). 

Case VI 

In the sixth case, we suppress the delta-function component {Z = 0) and we also impose 
the condition S{k = 0) = 0. We find that 0* = 0.3535 (see Tables 1 and 2). 



4 Concluding Remarks 

For the family of pair correlation functions specified by relation (8), the optimal packing 
fraction is characterized by unusual structural features such as substantially large density 
fluctuations in the infinite-wavelength limit, vanishing or nearly vanishing fluctuations at 
several finite wavelength values, and an interparticle radial distance (r ^ 1.4) at which particle 
centers are prohibited. Nonetheless, the maximal packing fraction and coordination number 
(0* = 0.627 and Z = 5.8) are consistent with values for dense random packings generated 
experimentally and computationally. Clearly, however, the properties (iv) and (v) ("split 
second peak" and the discontinuous drop at r = 2) that are characteristic of typical random 
packings are absent in the optimal solution. In future work, one may want to consider other 
families of functions that are not as smooth as (8) away from r = 1, e.g., piecewise continuous 
functions. Such extensions would quantify the significance of properties (iv) and (v) in raising 
0* above the optimal value of 0.627, while presumably driving S{k = 0) downward toward 
zero. 

Assuming that the optimal packing can actually be realized, there remain many open 
questions. Do the spheres form a contacting percolating network? Given the high density that 
is achieved, we suspect that the answer is in the affirmative. If so, what is the geometry of the 
contact set? Are rattlers present in the optimal solution? Is the packing locally, collectively, 
or strictly jammed?0 The answers to all of these questions can greatly be facilitated if we 
could determine whether there are packings that achieve the optimal solution. This can be 
accomplished using stochastic reconstruction techniques that enable one to obtain realizations 
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of sphere packings that have a targeted pair correlation function or structure factor. M We will 
attempt such a reconstruction in a future study. Another interesting extension of the present 
work is to generalize the family of pair correlation functions to the case of spheres of different 
sizes and to determine the maximal packing fraction. 
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Appendix: ti- dimensional Generalization of Case IV 

In this Appendix, we obtain an exact expression for the optimal value of (p* for the d- 
dimensional generalization of Case IV. A consequence of this result is a lower bound on 0* for 
random sphere packings in d dimensions, which we compare to a well-known lower bound for 
regular lattice packings. 



Consider the evaluation of the structure factor S{k) of relation (|T^ in d dimensions but 
without the short-ranged (damped-oscillating) contribution, i.e., 

S{k) = l + p[Gi{k) + Gnik)] (23) 

where Gi{k) and Gjj{k) are the d-dimensional Fourier transforms!! 

G,{k) = (2vr)t Hr'-' {,,(r) - 1} dr, (24) 

Gn{k) = {2.)i r r^-^ {gn{r) - 1} ^''''l^^fl} dr^ (25) 

g,(r) = U(t - 1). (26) 



and 
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The quantity Si{r) is the surface area of a d-dimensional sphere of radius r. 



It immediately follows that 



pGi{k) = YYXdE Jd/2{k), 



where 



d/2 



-{i/2r 



(29) 



(30) 



r(i + rf/2) 

is the (i-dimensional packing fraction. If this were the only contribution to the structure factor 
then the non-negativity condition S{k) > implies 



which agrees with the result given in Ref. 4. It easily follows that 



pGnik) 



d{kY/^-^ 



J, 



d/2-l 



Substitution of ( pOD and (|3^) into ( ^31) yields 



S{k) = 1 + {Z- 2U) + 



k^ + 0{k' 



1 + d/2 2d 

The last term changes sign if Z increases past 2'^(j)d/{d + 2). At this crossover point, 

2d+i 



Sik) 



+ 0{k' 



d + 2 

Since the minimum occurs at A; = 0, then we have the exact results 



Thus, we have the lower bound 



d + 2 



')* > 



2 



d + 2 



(31) 



(32) 



(33) 



(34) 



(35) 



(36) 



for random sphere packings in d dimensions because the addition of the short-range contri- 
bution (which we have neglected) would result in a generally larger value of (p*. In obtaining 
lower bound (|36| ) , we have assumed there are no further sufficiency conditions beyond (j^) and 
(I)- 

In very high dimensions {d ~ 1000), the densest known packings are non-regular lattices.ii 
Thus, it is of interest to compare the lower bound (|36| ) to the Minkowski-Hlawka theorem,!! 
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which gives a lower bound on the maximum packing fraction for d-dimensional regular lattices 
of identical spheres: 

(Pmax > (37) 

where C,{d) is the Riemann zeta function. The lower bound (|36|) is larger than the lower bound 
(p^) for (i > 3. Indeed, the difference between these bounds grows with increasing d. 
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Table 1: Terminal packing fractions (f)* for six different cases in which the step-function 
contribution gj to g in (|^) is always included. 



Case I 


911^0 


9III 7^ 


S{k) > 


0* = 0.627 


Case II 


911^0 


9III 


S{k) = 


0* = 0.46 


Case III 


911^0 


9III = 


S{k) > 


0* = 0.41 


Case IV 


911^0 


9III = 


S{k) = 


(p* = 0.3125 


Case V 


9ii = 


9III 7^ 


S{k) > 


(p* = 0.375 


Case VI 


9ii = 


9III 7^ 


S{k) = 


0* = 0.3535 



Table 2: Values of the parameters for the cases in Table 1. 



Case I 


A = 2.733 


B = 0.510 


C = 7.471 


D = 0.627 


Z = 5.80 


Case II 


A = 1.15 


B = 0.510 


C = 5.90 


L) = 1.66 


Z = 2.3964 


Case III 


A = 


B = 


C = 


D = 


Z = 3.1504 


Case IV 


A = 


B = 


C = 


D = 


Z = 1.5 


Case V 


A = 4.8 


B = 1.2 


C = 5.90 


D = 0.90 


Z = 


Case VI 


A = 3.9 


B = 0.9 


C = 5.70 


D = 0.90 


z = o 
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Figure Captions 



Figure 1: The pair correlation function g{r) vs. r as obtained by averaging over ten con- 
figurations at a packing fraction = 0.64. The "binned" peak value of g{l) (not shown) is 
approximately 24. 

Figure 2: Structure factor (a) and pair correlation function (b) for Case I. Note the appearance 
of a vertical line at contact in (b) indicating a delta-function contribution there. 

Figure 3: Structure factor (a) and pair correlation function (b) for Case II. Note the appear- 
ance of a vertical line at contact in (b) indicating a delta-function contribution there. 
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Figure 1 
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Figure 2a 
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A=2.733, B=0.510 
C=7.471, D=0.722 
Z=5.8, (|)=0.627 
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Figure 2b 
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A=1.15, B=0.510 
C=5.90, D=1.66 
Z=2.3964, ([)=0.46 
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Figure 3a 
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Figure 3b 
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